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ABSTRACT 

We present a numerical study of the cosmic density vs. velocity divergence relation 
(DVDR) in the mildly non-linear regime. We approximate the dark matter as a non- 
relativistic pressureless fluid, and solve its equations of motion on a grid fixed in 
comoving coordinates. Unlike N-body schemes, this method yields directly the volume- 
averaged velocity field. The results of our simulations are compared with the predic- 
tions of the third-order perturbation theory (3PT) for the DVDR. We investigate both 
the mean 'forward' relation (density in terms of velocity divergence) and the mean 'in- 
verse' relation (velocity divergence in terms of density), with emphasis on the latter. 
On scales larger than about 20 megaparsecs, our code recovers the predictions of 3PT 
remarkably well, significantly better than recent N-body simulations. On scales of a 
few megaparsecs, the DVDR predicted by 3PT differs slightly from the simulated one. 
In particular, approximating the inverse DVDR by a third-order polynomial turns out 
to be a poor fit. We propose a simple analytical description of the inverse relation, 
which works well for mildly non-linear scales. 

Key words: cosmology: theory - cosmology: dark matter ~ large-scale structure of 
the Universe - methods: numerical 



1 INTRODUCTION 

It is now widely believed that the large-scale structure 
formed by the growth of small inhomogeneities present in 
the early Universe. In this scenario, commonly referred to 
as the gravitational instability (GI) paradigm, cosmic density 
and velocity fields are tightly coupled, and the relation be- 
tween them involves the cosmological parameter fi. In the 
linear regime, i.e. for the r.m.s. density fluctuations much 
smaller than unity, the density-velocity divergence relation 
(DVDR) reduces to 



Blanton et al. 1998; Dekel & Lahav 1998) and observational 
(e.g., Davis & Geller 1976; Dressier 1980; Giovanelli, Haynes 
& Chincarini 1986; Santiago & Strauss 1992; Loveday et 
al. 1996; Hermit et al. 1996; Guzzo et al. 1997; Giavalisco et 
al. 1998; Tegmark & Bromley 1998) evidence that galaxies 
are biased tracers of the matter distribution. As a result, 
the comparisons between the fields in question within linear 
theory cannot yield an estimate of fl itself. What is actually 



measured is the quantity /3 
bias parameter. 



i} /b, where b is the linear 



5{r) = -/"'(n,A)V-v(r). 



(1) 



/(n,A)~n°-V 



(2) 



Here, S is the mass density fiuctuation field, v is the peculiar 
velocity field, distances are expressed in units of km s~^, and 

70 V 2 J 

(Lahav et al. 1991). The factor / depends mainly on Q and 
only weakly on the cosmological constant A (provided that 
A is in the range allowed by observations). The comparisons 
between density and velocity fields are a useful test of the 
GI hypothesis. In principle, they may also be used as a tool 
to measure fl (Dekel et al. 1993). 

However, there is both theoretical (e.g., Kaiser 1984 
Davis et al. 1985; Bardeen et al. 1986; Dekel & Silk 1986 
Cen & Ostriker 1992; Kauffmann, Nusser & Steinmetz 1997 



The current state of estimates of l3 is confused. The so- 
called velocity-velocity comparisons generally result in low 
values of (3 (~ 0.5: Roth 1994; Schlegel 1995; Schaya, Peebles 
& TuUy 1995; Davis, Nusser & WUlick 1996; da Costa et 
al. 1997; Riess et al. 1997; Willick et al. 1997; Willick & 
Strauss 1998), while density-density comparisons yield high 
values (~ 1.0: Dekel et al. 1993; Hudson et al. 1995; Sigad 
et al. 1998) . In velocity-velocity comparisons, galaxy density 
field is used to predict the associated peculiar velocity field, 
which in turn is compared to the observed peculiar velocities 
of a sample of galaxies with measured redshift-independent 
distances. In density-density comparisons, velocity data are 
used to reconstruct the underlying mass density field, in 
order to compare it with an observed galaxy density field. 
A number of possible explanations of the divergence in the 
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estimated values of (5 has been proposed (see, e.g., Sigad et 
al. 1998). One of them are non-linear effects. 

The density fluctuations obtained from current redshift 
surveys (e.g. Fisher et al. 1994) and from the potent (Dekel 
et al. 1998) reconstruction of the mass density field slightly 
exceed the regime of applicability of linear theory. For ex- 
ample, the density contrast in regions like the Great Attrac- 
tor or Perseus-Pisces is around unity even when smoothed 
over scales of 1200 km s~^, currently employed in density- 
density comparisons (Sigad et al. 1998). In velocity-velocity 
comparisons, the fields in question are generally smoothed 
over smaller scales than in density-density ones. Aston- 
ishingly, while in current density-density comparisons the 
non-linear corrections to the linear density-velocity relation, 
equation (^j), are accounted for, in velocity-velocity com- 
parisons they are not. The only exception is an attempt by 
Willick et al. (1997) to model the DVDR by a second-order 
formulaFJ To their surprise, the maximum-likelihood fit of 
the predicted to the observed peculiar velocities was for zero 
amplitude of the second-order corrective term. However, the 
smoothing scale they used was 3h~^Mpc. At such a small 
scale, the variance of the density field is already in excess 
of unity and, as we will show later, neither the linear nor 
the second-order formula is a good description of the actual 
DVDR. 

The purpose of this paper is to propose a simple and ac- 
curate description of the DVDR at mildly non-lineair] scales, 
which would be easy to implement in current velocity- 
velocity comparisons. To date, there have been several at- 
tempts to construct a mildly non-linear extension of re- 
lation (nl). They were either based on various analytical 
approximations to non-linear dynamics (Bernardeau 1992; 
Catelan et al. 1995; Chodorowski 1997; Chodorowski & 
Lokas 1997, hereafter CL97; Chodorowski et al. 1998, here- 
after CLPN), or N-body simulations (Mancinelli et al. 1994; 
Ganon et al., in preparation), or both (Nusser et al. 1991; 
Gramann 1993; Mancinelli & Yahil 1995). So far, the most 
comprehensive description of the mildly non-linear DVDR 
has been recently done by Bernardeau et al. (1999; hereafter 
B99). Our work is an extension and improvement of B99 in 
several ways: 

• In B99 the analysis was for technical reasons performed 
solely for fields smoothed with a top-hat filter. Here we also 
analyze fields smoothed with a Gaussian filter, which is now 
commonly applied to observational data. 

• The fully non-linear formula proposed by B99 expresses 
density in terms of the velocity divergence (the so-called 'for- 
ward' relation). However, in velocity-velocity comparisons 
one needs a formula for the velocity (divergence) expressed 
as a function of the density ( the so-called 'inverse' relation). 
Due to the scatter in the DVR, the latter is not given by a 



straightforward inversion of the former. We obtain such an 
'inverse' formula here. 

• Our 'inverse' formula is much simpler compared to the 
'forward' one of B99, but equally accurate, as detailed com- 
parisons with numerical simulations show. Unlike the second 
order formula used by Willick et al. (1997), it works well for 
smoothing scales down to a few megaparsecs. 

• Instead of performing N-body simulations, we model 
cold dark matter as a pressureless cosmic fluid. We solve 
non-linear equations for its evolution on a grid fixed in co- 
moving coordinates. This approach is advantageous over the 
standard N-body one for studying the evolution of the ve- 
locity field in the mildly non-linear regime. The reasons are 
outlined below. 

Both in N-body simulations and in our code, the final 
velocity field is known at a discrete set of points. In the case 
of an N-body simulation this set is particles' positions, in 
our case it is the grid. Due to clustering, the N-body veloc- 
ity field is sampled very non-uniformly, while the sampling 
of our velocity field is perfectly uniform. Smoothing of a 
non-uniformly sampled velocity field leads to the so-called 
'sampling gradient bias' (Dekel, Bertschinger & Faber 1990). 
In N-body simulations, the sampling rate of the velocity field 
is proportional to the number density of particles in a given 
region. The averaging of the field within a smoothing win- 
dow is therefore not volume- but mass-weighted, resulting in 
a special type of bias mentioned above. To circumvent this 
problem, elaborate 'tessalation' algorithms for the velocity 
field have been proposed (Bernardeau & van de Weygaert 
1996). However, they work only for a top-hat filter. Another 
problem is that N-body simulations provide very little infor- 
mation on the velocity field in voids, simply because there 
are very few velocity tracers there. 

Due to uniform sampling, our simulations yield directly 
volume-weighted values of velocity, for any type of smooth- 
ing. Moreover, we probe the velocity field in the voids as 
finely as in dense regions. As a result, at a very low numeri- 
cal cost it was possible to have the velocity field sampled at a 
comparable number of points to that of B99 (64^ compared 
to 50^), and still of significantly better quality, as shown 
below. 

The paper is organized as follows: In Section Gi we dis- 
cuss the theoretical aspects of the DVDR. Then, in Sec- 
tion M, we present our simulations, describing the algorithm 
in Section 3T^ and the cosmological model investigated in 
Section 3.2. In Section W we investigate the mean forward 
density-velocity relation and we demonstrate that it is w ell 
described by a third-order polynomial. In Section 



5.1 



show that the polynomial formula is a poor approximation 
of the inverse relati on, and we propose an alternative de- 
scription in Section 5.2, We summarize our results in Sec- 
tion |6[ 



* Strictly speaking, they proposed a fully non-linear formula, but 
in the process of actual comparison they truncated it at second 
order terms. 

1 We define mildly non-linear scales as these at which the r.m.s. 
density fluctuation is a significant fraction of, but still smaller 
than, unity. Then the mildly non-linear scales in the Universe 
are about or greater than 8 h~^Mpc for top-hat smoothing, and 
roughly twice smaller for Gaussian smoothing. 



2 THEORETICAL FRAMEWORK 

Due to the Kelvin circulation theorem, the cosmic veloc- 
ity field remains irrotational before shell-crossings. It can 
therefore be described by a single scalar function, which we 
choose here to be the velocity divergence, V ■ v (through- 
out the paper the derivative is taken in velocity units, i.e. 
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8 Mdc top hat filter 



Figure 1. Joint probability distribution function of density and 
velocity divergence, P{5, V ■ v), from a 64^ simulation. Data are 
convolved with a 81i~^Mpc top-hat filter. Contours have intervals 
of 0.5 in logj^g(P). Solid line represents the linear relation. 

H = 1). The linear relation (eq. Il|) between the density con- 
trast and the velocity divergence at a given point holds only 
on scales large enough so that the density fluctuations are 
much smaller compared to unity. On smaller scales, non- 
linear effects modify the relation in a number of ways. For 
full discussion of the density versus velocity divergence re- 
lation in the mildly non-linear regime the reader is referred 
to B99. 

In brief, qualitative features of the relation can be out- 
lined as follows: 

• It is non- linear. 

• It is also non-local, which implies that it is locally non- 
deterministic, i.e. it has a scatter in S for a given V • v and 
vice versa. 

• Since the scatter originates exclusively from higher 
(than linear) order terms, it is small in the mildly non-linear 
regime. Therefore, the most probable values of S and V ■ v 
form an elongated region in the (5, V ■ v) plane. 

Figure nl is a typical plot of the values of S and V ■ v 
obtained in our simulations. The fields are smoothed with 
a top-hat filter with the smoothing radius of 8 h~^Mpc. 
For such a smoothing scale, the r.m.s. fiuctuation of the 
density field, as, in our simulations is ~ 0.9, so the fields are 
close to leave the regime of mild non-linearities. However, in 
Figure hi one can still observe an obvious correlation between 
the density and velocity divergence. 

Analytical calculations predict (Bernardeau 1992; Gra- 
mann 1993; Catelan et al. 1995; MancineUi & Yahil 1995; 
Chodorowski 1997; CL97; CLPN) and N-body numerical 
simulations confirm (MancineUi et al. 1994; B99; Ganon et 
al. 1999) that the mildly non-linear DVDR depends on fi 
and A in a very simple way. Specifically, if we define the 
scaled velocity divergence. 



the relation between the density and the scaled divergence 
(for simplicity it will also be referred to as DVDR) will be 
practically Q- and A- independent. Since the relation has a 
scatter, the full information about the DVDR is contained in 
the joint probability distribution function (PDF) for 5 and d. 
Such a joint PDF has been constructed by B99. However, the 
scatter is small compared to random errors in the observed 
density and velocity fields (CLPN). Therefore, of most in- 
terest for practical applications are the mean relations: the 
mean density for the given velocity divergence, (S^g) (the 
'forward' relation), and vice-versa, (^i^) (the 'inverse' rela- 
tion) £| The forward relation is relevant for density-density 
comparisons; the inverse relation is relevant for velocity- 
velocity comparisons. Since velocity-velocity comparisons 
employ smaller smoothing lengths, non- linear effects are 
more important there than in density-density comparisons. 
That is why in this paper we shall concentrate on finding a 
simple, and simultaneously robust, description of the inverse 
relation for Gaussian smoothing of the fields. 

Though one might formally derive the inverse relation 
from the joint PDF constructed by B99, it would be inap- 
propriate for a number of reasons. Firstly, while the forward 
relation can be derived from this PDF in an analytic form, 
the inverse one can only be computed numerically. Secondly, 
the joint PDF was constructed by B99 for top-hat smoothed 
fields and it is expected to depend quantitatively on the type 
of smoothing. Finally, with our fiuid code we hope to trace 
the actual DVDR more accurately. 

The mildly non-linear regime is the one in which per- 
turbation theory can be applied. In particular, the mean 
relations are a priori accessible to analytical perturbative 
calculations. CL97 derived the forward DVDR up to third- 
order terms, accounting for the smoothing of the density and 
velocity fields. The mean density contrast given the scaled 
velocity divergence is a third-order polynomial in the diver- 
gence, 

{S^g) = ao + aie + a2e^ + aae\ (4) 

where ao = —02 erg and ag is the variance of the scaled veloc- 
ity divergence field, {d'^). The coefficients, a^, appearing in 
the above expansion were explicitly calculated by CL97 for 
Gaussian smoothing and by B99 for top-hat smoothing. As 
explained above, they depend extremely weakly on Q and A. 
CLPN derived the inverse relation up to third-order terms. 



5) = ro + riS + r25 + rgS . 



(5) 



The coefficients Vi were calculated by CLPN for Gaussian 
smoothing and by B99 for top-hat smoothing. 

Contributions to the DVDR from orders higher than 
third are known in only one special case, of unsmoothed 
fields with vanishing variance. Bernardeau (1992) derived 
for this case the following formula: 



,) = |[(l+5)2/3„l] 



(6) 



-ri(n,A)v-v, 



(3) 



The above expression is strictly valid only for ct^ — > and 
n ^ 0, but since the il-dependence of the (scaled) DVDR 



■f Due to the scatter, the inverse relation is not given by a 
straightforward inversion of the forward one. 
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Figure 2. P{9), the probability distribution of the scaled ve- 
locity divergence for 8h~^Mpc top-hat smoothing. Open squares 
are combined results from our six runs, with error-bars shown. 
Solid line represents formula (12) of Bernardeau (1994), with the 
variance of 9 taken to be the average from our simulations. 



is extremely weak, it remains a good approximation also for 
other values of fl. 

Equations (O) and m are two difTerent approximations 
to the inverse relation. As already stated, equation (H) ac- 
counts for smoothing and for finite variances of the fields. 
Equation (a) does not,p|but instead it includes contributions 
from all orders. We can therefore expect the two equations to 
carry complementary information about the actual relation. 
Our procedure of finding a simple and accurate description 
of the inverse relation will consist of two steps. Firstly, we 
will check on which scales the third-order expression (Mj is 
a good description of the relation, and at which scales it al- 
ready fails. Then, guided by our numerical results, and by 
equation (h), we will look for a formula for the inverse re- 
lation, which would be accurate in the whole range of the 
mildly non-linear scales. 



3 THE SIMULATIONS 

3.1 The code 

As stated in Section hi, despite their simplicity and numer- 
ous advantages, N-body simulations of cosmic velocity fields 
have several drawbacks. 

To cope with them, we have performed our simulations 
using CPPA (Cosmological Pressureless Parabolic Advec- 
tion): our original Eulerian, uniform-grid based code for self- 
gravitating pressureless fluid evolution in an expanding Uni- 
verse. The main ideas of the algorithm are similar to those 



3 B99 argued however that this result should remain valid for 
top-hat smoothed fields with vanishing variance. 



of Peebles (1987), with several modifications. An early ver- 
sion of the code is described in Kudlicki, Plewa & Rozyczka 
(1996); later improvements to the code and test results will 
be described in detail in Kudlicki et al. (in preparation). 
CPPA employs a three-dimensional Cartesian grid fixed in 
dimensionless comoving coordinates (see Gnedin 1995). The 
Poisson equation is solved by a standard FFT-based routine, 
working on the same grid as the Euler solver. Advection of 
mass and momenta is done using the piecewise-parabolic 
scheme, as described by Colella & Woodward (1984). The 
advection step consists of a series of sweeps along the main 
axes of the computational domain. Unlike Peebles (1987), 
we use a variable timestep, according to the CFL condi- 
tion. To allow for shell crossing on small scales, and inhibit 
the unrealistically high densities at cluster centres, we ar- 
tificially interchange fiuxes across local directional maxima 
of the density field. During a sweep, if a local directional 
maximum of density is encountered, and there is matter 
falling onto it from both sides, then the fiuxes of density 
and momentum calculated at the left and right interface of 
the maximum density cell are interchanged. Performed in all 
directions, it successfully inhibits the non-physical transfer 
of power into the smallest scales, while conserving the total 
momentum. We are aware that our code does not reproduce 
small-scale structures properly, however for the present pur- 
pose - involving window functions larger than 3 h~^Mpc 
(Gaussian) and 6 h^^Mpc (Top Hat) - it is an efficient and 
satisfactory tool. 

3.2 Selection of the parameters and the models 

Since the relation between density and scaled velocity di- 
vergence depends very weakly on the background cosmo- 
logical model (see the previous section), we were free to 
choose the convenient and well-tested Einstein-de Sitter 
model (n = 1,A = 0). 

To estimate random errors, we have performed six real- 
izations of this cosmological model, each of them with dif- 
ferent random phases of the initial density field. 

We aimed at investigating the statistics of density and 
velocity fields on both intermediate (several megaparsecs) 
and large (up to 60h~^Mpc) scales, so, in order to sup- 
press the effects of finite simulation volume size and improve 
the statistics we decided to make the simulation box signifi- 
cantly larger than the largest filter used. A reasonable solu- 
tion turned out to be a (200 h~^Mpc)'' cube with standard 
periodic boundary conditions. For 64^ grid cells, the spa- 
tial resolution is 3.125 h~^Mpc, sufficient for our purposes. 
As a test we performed a simulation with 128'' grid cells 
for a (200h~^Mpc)^ cube and the results remained in good 
agreement with those obtained with the coarser grid. 

Most of the analytical calculations in the discussed 
regime have been done for scale-free, or power-law, power 
spectra, P{k) ex k" . These spectra may seem artificial, 
none the less they are believed to approximate the real 
power spectrum at least piecewise over significant ranges 
of wavelengths. In particular, in the range of scales ^ 1- 
20h~^Mpc, the observed power spectrum is well approxi- 
mated by a power law (e.g., Sutherland et al. 1999; Freudling 
et al. 1999). That is why we have decided to use a power law 
power spectrum in our simulations. For convenience, we have 
picked such normalization of the initial density contrast. 
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Figure 3. The coefficients 01,02,03 from six simulations (curves and error-bars) and their third-order perturbation theory predictions 
(solid lines). Left panel: Gaussian filter, right panel: Top-Hat filter. 
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Figure 4. The coefficients r'i,r2,r3 from six simulations (curves and error- bars) and their third-order perturbation theory predictions 
(solid lines). Left panel: Gaussian filter, right panel: Top-Hat filter. 



that the amphtude of the hnear growing mode measured 
with a Sh^^Mpc top-liat filter at present epoch is unity, 
(1 -I- 2initiai)o"8, initial ~ 1- Siuce we have chosen a scale-free 
power spectrum, the results may be simply rescaled to other 
normalizations. 

Observations suggest that on mildly non-linear scales 
the effective spectral index n lies between —1 and —1.5 
(Baugh & Efstathiou 1993, 1994; Fisher et al. 1993; Feld- 
man et al. 1994; Park et al. 1994; Lin et al. 1996; Suther- 
land et al. 1999). For our simulations we have chosen the 
value of n = — 1, because one of our code tests was to com- 
pare the PDF of the velocity divergence with the analytical 
formula of Bernardeau (1994). This comparison was essen- 
tial to demonstrate the advantage of our code over N-body 
simulations in velocity field studies. 

N-body schemes yield mass- weighted velocity fields, re- 



sulting in spurious velocity gradients. These gradients man- 
ifest themselves as spurious tails in the PDF of 6. Tesse- 
lation techniques, invented by Bernardeau & van de Wey- 
gaert (1996) to overcome this problem, are very CPU-time- 
consuming, and the results have much lower resolution than 
the simulations themselves. (Commonly 50"^ compared to 
128'^: Bernardeau & van de Weygaert 1996; Bernardeau et 
al. 1997; B99). In contrast, our code yields directly a volume- 
weighted velocity field. Figure presents the comparison of 
the PDF of 6 as recovered from our numerical data with the 
analytical formula of Bernardeau (1994). The velocity field 
is smoothed with a top-hat filter of the radius of 8 h~^Mpc. 
Note that the PDF obtained from the simulations has no 
spurious tails, whatsoever. The definition of the scaled di- 
vergence is such that negative 6 corresponds to positive di- 
vergence, i.e., to the expansion of voids. Note how well the 
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negative tail of our PDF traces the analytical prediction, on 
over three decades of the probability value. (Compare Fig. 6 
of Bernardeau 1994.) In the extreme part of the positive 
tail, the analytical PDF slightly overestimates the measured 
one. This is to be expected, since the formula of Bernardeau 
(1994) is only an approximate fit to the actual PDF, overes- 
timating the value of skewness of the distribution (2 instead 
of 1.7 based on PT). Indeed, in our simulations the skewness 
measured with an 8 h~^Mpc top-hat filter has the value of 
1.77 ±0.11. 



4 THE FORWARD RELATION 

With our models, we have tested the polynomial approxi- 
mation of the mean forward DVDR (Eq. kh. 

We compare the coefficients ai , 02 and as computed 
from our simulations to the corresponding third-order PT 
values in Figure y. These coefficients have been computed 
from the simulations using a standard four-parameteiPI least 
square fit on data from each of the 64'^ grid cells (after 
smoothing with the filters required) . We have also performed 
a fit for all the coefficients ao through 05 , in this case a4 and 
as were consistent with zero, and the values of ao ■ ■ -0,3 did 
not change remarkably between these two fits. With the ex- 
ception of the most highly non-linear smoothing scales, the 
values of a„ very weakly depend on the filter size, and are 
in good agreement with the perturbation theory predictions 
of CL97. 



5 THE INVERSE RELATION 

5.1 Polynomial parameterization 

As we have shown above, the mean forward relation can 
be described with sufficient accuracy by the polynomial for- 
mula (W) on the relevant scales. In this section we test the 
polynomial approximation (H) to the inverse relation. 

In Fig. W we present our numerical estimates of the pa- 
rameters ri, r2 and ra. On scales below ~ 8 Mpc for Gaus- 
sian and ~ 15 Mpc for top-hat filtering, their dependence on 
the filter size is very strong, especially for r2, which is the 
first, and therefore the most essential parameter describing 
the non-linearity of the DVDR. This makes the polynomial 
formula inconvenient for application to observational data. 
Even using the value of r2 predicted for a particular filter 
size will not help much since all sizes scale with the present- 
epoch density contrast, ag, and the Hubble constant, h, and 
neither of these parameters is known accurately yet. 

We have also performed a fit for the parameters ro 
through rs of the fifth-order polynomial in S. The values 
of r4 and rs turn out to be inconsistent with zero at more 
than 2a level for filter sizes smaller than 20 h~^Mpc, and, 
which is still more important, their addition to the fit signif- 
icantly infiuences ri, r2 and rs (see Figure m. Moreover, a 
third-order polynomial that fits the distribution in the large- 
scale regime (for small \5\) obviously has a non-monotonic 



We have found that the value of ao is perfectly consistent with 
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Figure 5. The value of r2 from 4-parameter fit (thin solid curve) 
and 6— parameter fit (dashed curve). Heavy solid line represents 
the third-order perturbation theory value. 
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Figure 6. The a parameter of the (Oig) relation. Dotted curve 
and error-bars: fit with second-order approximation to e; solid: a 
from fit of both a and e. 



derivative, which is inconsistent with the properties of the 
actual distribution. Figure \u clearly demonstrates that the 
first derivative of the relation is positive and monotonically 
decreasing. The third-order polynomial fit to the simulated 
data is drawn with a long-dashed line in Figure m. One can 
observe that this function has an inffection point well within 
the range of 5 and 9 occupied by numerical data, not seen in 
the simulated relation. Another substantial disadvantage of 
high-order polynomial fits is the high number of parameters, 
each of them depending on the smoothing scale. Finally, the 
Tn parameters are strongly correlated, so possibly another 
formula, dependent on fewer parameters and not a polyno- 
mial, will provide a better description. 



— a2(T^, and that a 3-parameter fit with the ao 
straint gives the same values for ai , a2 and as . 



-a2crfj con- 
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Figure 7. Joint PDF of 5 and 0, P{S,6): contours liave intervals of 0.5 in log]^Q(P); combined data from six simulations are plotted. 
Solid curve represents our formula M) with a and e fitted independently, dotted curve - a one-parameter fit according to formula (VV) 
with the e offset accurate to the second order (eq. Sh . Long-dashed curve represents the third-order polynomial fit (H) , and short-dashed 
— the non-linear formula of Willick et al. (1997). Window in the upper-left corner enlarges the void region, here for clarity contours have 
intervals of 1.0 in log]^Q(P). 



5.2 A robust non-linear formula 

In search for a better formula for the (^i^) relation one needs 
to take into account: (a) monotonicity of the fit and its 
derivative, (b) agreement with the polynomial description 
for large filter radii, (c) proper asymptotic behavior in the 
voids, and (d) the mass conservation law, i.e. the formula 
should yield (9) = 0. Of the above, (a), (b) and (c) are sat- 
isfied by the asymptotic formula of Bernardeau (1992), (see 
eq. o of this paper), which was derived in the limit of zero 
density dispersion, Of, —> 0, and which does not account for 
smoothing effects. Filtering of the data substantially affects 
the higher order moments of the distribution of 5 and 8, 
and therefore it is expected to cause a change in the shape 
of the DVDR. In order to make equation m applicable to 
filtered data (and such are all the data obtained from galaxy 



catalogs) we have made an educated guess, replacing the ex- 
ponent of 2/3 in the formula with a free parameter, l/Q.fll 
This change does not affect the shape of the fit in the void 
wing of the plot, as long as a is not very much different 
from |. To satisfy (d), i.e. to keep the average 6 equal zero, 
we add a constant, depending only on a and scaled by the 
density dispersion. Our final formula has the form: 



= q[(1 + <5)'/°-1] 



(7) 



where the constant e can be approximated as: 



II This follows the idea of B99 in a sense, but will result in a much 
simpler formula. 
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2q 



1 2 



(8) 



Formula (n) is accurate to the second order, we have tested 
its relevance by fitting both a and e as independent parame- 
ters. The values of a obtained in these two fits are consistent 
for filter radii larger than about 6h~^Mpc (see Figure |6|). 

We plot the joint PDF of 5 and 6 combined from all our 
six simulations in Figure M. Willick & Strauss (1998) per- 
formed their VELMOD analysis of peculiar velocities using 
IRAS galaxy density field convolved with either 3h~^Mpc or 
5h~^Mpc Gaussian filter. They reported very similar results 
for both smoothing scales. Therefore, we chose to plot the 
distribution for data fitered with 4h~^Mpc Gaussian kernel. 

Formula nm is drawn in this figure with heavy solid 
line (two-parameter fit), and dotted line (one-parameter fit, 
offset accurate to the second order, formula H). The long- 
and short-dashed lines represent respectively the third-order 
polynomial fit and a formula by Willick et al. (1997): 



5) = 



(l + aV|)^ + acr| 
1 + a5 



(9) 



in which a is a free parameter. The polynomial fit is poor 
not only for very large 5, but it also overestimates 9 for 
1 < 5 < 5. The formula of Willick et al. follows the PDF 
as closely as ours for —0.6 < 5 < 4, but departs from the 
simulated distribution at the extreme values of 5. We have 
also estimated a from our simulations, obtaining the value 
of 0.14 at the 4h~^Mpc scale. At large scales our estimate 
of a is rather 0.24 than 0.28 reported by Willick et al. 

Our fit for the one-parameter a-formula slightly over- 
estimates 9 in the high-(5 tail. With the two-parameter de- 
scription used instead, our formula fits the data very well in 
the entire range of (5, at a slightly greater a. 

The weak dependence of ct on the smoothing scale is well 
visible in Figure o (in the two-parameter description, a ~ 
1.9 ±0.1). As stated earlier, the formula (fTI) as a description 
of the inverse relation was our 'educated' guess. Expand- 
ing it for large smoothing scales, and comparing it with the 
polynomial description we find that ct ~ oils = 1/(1 + 2r2). 



From F gure 4 we obtain ctLS — 1.95, which remains in good 
agreement witti direct tit. i^br very small smoottimg scales, 
simple considerations of energy conservation in the model of 
spherical collapse yield a = 2. We expect a to change very 
weakly between weakly and highly non-linear scales, which 
is indeed observed. 

Our formula is also much simpler than formula (18) of 
B99 for the forward relation. The reason that the formula 
of B99 is complex is twofold. Firstly, B99 aimed at mod- 
elling the weak ^-dependence of the relation. However, the 
f2-dependence turned out to be so weak that it was practi- 
cally unnecessary to account for it. Secondly, they rigorously 
applied the constraint, coming from the maximal expansion 
of voids, that 9 = —3/2 for 5 = — 1. Although we have not 
required it explicitely, our formula satisfies this 'voids' con- 
straint very well (see Fig. 7). 



6 SUMMARY AND CONCLUSIONS 

We have tested several parameterizations of the density vs. 
velocity divergence relation on weakly and mildly non-linear 
scales. We confirm that the polynomial formula provides a 



good description for the forward relation (the density con- 
trast as a function of the velocity divergence) . 

On the other hand, the inverse relation is not well de- 
scribed by the polynomial expansion, which does not con- 
verge fast enough. Also, on mildly non- linear scales the pa- 
rameters of the expansion strongly depend on the smoothing 
scale. The formula of Willick et al. (1997) is better than the 
polynomial description, but it is not free from drawbacks, 
either. Firstly, for large densities it has a horizontal asymp- 
tote, not observed in the simulated distribution. As a result, 
in the high-density tail it underestimates the actual relation. 
Secondly, like the polynomial description, it incorrectly de- 
scribes the low-density tail (i.e., the relation for voids). 

Our formula (7) is free from these disadvantages. It it 
also simpler (in its simplest form it depends on one parame- 
ter only) , which makes it easy to implement in the velocity- 
velocity comparisons. 
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